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ABSTRACT 


Signatur e 




THE TIME DEPENDENT FINITE ELEMENT MODELING OF THE ELECTROMAGNETIC 


FIELD IN ELECTRICAL MACHINES - METHODS AND APPLICATIONS 


Elias G. Strangas, Ph.D. 

University of Pittsburgh 

A method is developed to calculate the field quantities at a cross 
section of an electrical rotating machine, and to use these quantities 
to calculate the machine parameters and performance. The method includes 
the effects of eddy currents and saturation and calculates the inductances, 
torque and losses due to eddy currents. It is applicable to cases which 
cannot be described by a constant set of parameters, but rather where 
these parameters become functions of time and field solution. 

The method of time dependent finite elements is used, and a technique 
is developed to describe the airgap when the rotor is moving with respect 
to the stator. The solution techniques that can be applied to the resulting 
system of equations are investigated, and the method of the preconditioned 
conjugate gradient is described and applied. The underrelaxation method 
is applied to the calculations of the saturation levels in iron parts. 


iii 





The solution of the magnetic potential at the cross section of the 
machine is used to calculate inductances through flux linkages, and torque 
through the application of the Maxwell stress tensor at a surface in the 
airgap. The techniques developed are used to calculate the performance 
of chopper controlled DC series machines used in electric vehicles and 
to predict torques and currents during the asynchronous starting of 
synchronous salient pole motors. 
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1.0 INTRODUCTION 

Availability of modern, high speed digital computers with large 
storage capacity and speed of operations has had a dramatic impact on 
the methods of treating electrical machinery problems. Sophisticated 
numerical techniques, yielding accurate results, have replaced the 
approximate analytical approaches of the past. 

The modern techniques are based upon rigorous mathematical analysis 
methods; also they have stimulated the further development of the methods. 
A case in point is the problem of calculating the electromagnetic field 
of a machine. Initially only the field in small regions was calculated 
by rather crude means based on the fundamental magnetomotive force and 
flux waves. Then, the finite difference method emerged, providing a 
powerful tool with capability and accuracy exceeding that of the earlier 
methods. Following this, the finite-element approach to machine analy¬ 
sis was introduced, providing a still more capable analysis technique. 

In addition to the fact that numerical solutions are more elegant, 
fast, and accurate than analytical calculations or hand flux mapping, 
they are able to yield solutions to new problems which, because of their 
complexity, could not be treated with the old methods. A typical example 
discussed in this dissertation, is that of the performance of motors 
operated from a controller which produces a complex voltage waveform. 

An analytical treatment of such a problem requires many simplifications 
which introduce relatively large errors and which do not provide accu¬ 
rate knowledge of the internal conditions of the machine. The problem 
of starting a synchronous motor, although entirely different, has the 
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same characteristics as the operation of a machine from a controller, 
since the varying saturation and speed result in continuously changing 
operating conditions and parameters. 

1.1 Numerical Calculation of the Electromagnetic Field 

The first numerical techniques applied to the solution of the 

electromagnetic field in a machine were derived not from a mathematical 

theory which would provide an approximation to the exact solution of the 

field, but instead from the development of an equivalent circuit that 

used capacitors and resistors to represent conductivities and reluctances. 

The voltages calculated from this circuit were equivalent to the value 

of the magnetic vector potential in a cross section of the machine. 

Later, the mathematical method of finite differences was used to provide 

the systems of equations that would yield the field quantities. Erdelyi 
(1 2)* 

and Fuchs ’ presented a series of papers using this method for the 

calculation of the field distribution in DC machines, accounting for 

saturation in the iron regions. They later applied their technique to 

( 3 - 5 ) 

synchronous machines. This approach resulted in flux plots from 

which self and mutual inductances could be calculated as well as steady 

(6 71 

state, transient and subtransient reactances. Demerdash and Hamilton'' ’ 
demonstrated that by using an effective permeability which yields the 
same energy in the field as does the integral /BdH over one half of a 


^Parenthetical references placed superior to the line of text refer 
to the bibliography. 
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period, one may assume a sinusoidal flux density and magnetomotive force 

(O') 

and solve the problem in the frequency domain. Salon v J used the finite 

difference method to calculate eddy currents in generator rotors by 

solving for the field distribution in the time domain. He used the 

solution in the frequency domain to calculate voltages, currents and 

inductances and to solve a plethora of problems related to the operation 

of synchronous cylindrical rotor generators. 

( 9 ) 

In 1971, Chari introduced the finite element method to the sol¬ 
ution of electromagnetic fields. Although much more complicated to pro¬ 
gram than the finite difference method, it was proved superior for the 
modeling of boundaries and contours of materials and for giving detailed 
information about the field in specific areas of interest. Demerdash 
and Niehl^^ demonstrated in 1977 the superiority of this method in 

accounting for the saturation of iron. In 1975, Foggia, Sobonnadiere and 

( 11 ) 

Silvester presented a time dependent solution for the saturated eddy 
current case. A theoretical analysis of the time dependent finite ele¬ 
ment method is found in the thorough treatment of the subject by Douglas 
( 12 ) 

and Dupont , where the nonlinearity is taken into consideration with 
the use the of Crank-Nicolson-Galerkin approximation. 

1.2 The Calculation of Voltages, Inductances and 
Saturation of Iron Domains 

The solution of the electromagnetic field in an electrical machine 
requires knowledge of the current densities in the conductors. This 
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presents difficulties since what is known in most cases are the voltages 
at the terminals of the machine, rather than the currents. When the 
performance of the machine is to be determined, one has to resort to the 
use of self and mutual inductances. These are calculated by injecting 
currents into the conductors, solving the electromagnetic field and from 
that calculating the voltages induced in all of the windings. In DC 
machines this is quite straight forward. In AC machines it is simpler 
to use the "two Reaction" Theory and to calculate d- and q- axis para¬ 
meters for various current and saturation levels. This is a model 
developed from the Blondell-Park transformation and results in a system 
of differential equation, which can be used to predict the machine 
performance. 

The above approach contains certain unavoidable approximations. 
Firstly, the inductances are considered constant, whereas they actually 
depend upon the saturation level and thus are functions of the currents. 
Secondly, the two reaction theory presupposes linearity, so that super¬ 
position can be possible, which is not the case when saturation is 
present. Thirdly, the effect of eddy currents and skin effect on the 
inductances and the field solution is neglected, but usually this cannot 
be done without significant loss of accuracy. 

In the time domain solution, one can do away with inductances by 

incorporating in the system of equations resulting from the finite 

element method, the dependence of the current on the applied and induced 

( 8 ) 

voltages^ . This technique, however convenient, increases drastically 
both the storage requirements and the number of the operations needed 
for the solution of the resulting system of equations. 
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1.3 The Asynchronous Starting of Synchronous Motors 

The main consideration in the design of a synchronous motors is 
obviously the performance during synchronous operation, while the 
starting performance often receives less attention. A synchronous motor 
is usually started by connecting across the line or through an auto 
transformer. Only in the case of very large motors are they brought up 
to near synchronous speed by a smaller "pony" motor before the armature 
is energized. In the case of direct starting, the field excitation 
terminals are not energized, but are connected through a resistor of 
comparatively large ohmic value for the starting up period and the field 
is energized at near synchronous speed. The torque developed during the 
start-up is that of an asynchronous motor, but due to the dissymmetry of 
the rotor it contains an oscillation of twice the slip frequency (120-0 
Hz for a supply of 60 Hz). These torque oscillations can be in resonance 
with the connected load, resulting in severe stresses on the shaft and 
gear teeth used to connect motor and load. Also, the currents induced 
in the damper bars may increase the temperature to the point that the 
bars and pole faces are overheated and distorted. 

These problems have been recognized and addressed for a number of 

(13) 

years. As early as 1930, Linville wrote a classic paper in which he 
derived equivalent circuits for rotor and stator; resolved the starting 
currents and magnetomotive forces into forward and backward components 
and applied the Blondell-Park equations to calculate currents and torque. 

Since then, numerous studies have been conducted to cover both the 
mechanical and electrical aspects of asynchronous starting^^ . The 
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( 26 ) 

majority are along the lines of Linvill's work. Barret^ J recently 
(1980) presented a technique for predicting the starting performance of 
a motor with solid, salient poles, by calculating, in the time domain, 
the magnetic field using the finite element method in the airgap and an 
analytical one in the rotor iron. 

The technique of developing equivalent circuits for the rotor and 
stator, in the d- and q- axis, has been refined with the use of numeri¬ 
cal tools, and more lumped elements have been added; some, time and slip 
dependent to account for the space harmonics of the magnetic field and 
skin effect. Goodman'" J included higher current harmonics and corres¬ 
ponding inductances in the model of the machine for which he calculated 

(27-30) 

the torque during asynchronous starting. The work of Jovanovski 
strongly supported by experimental results, gives an insight in the cur¬ 
rent distribution in the damper winding, both in the case of squirrel 
cage and grill connections of the damper bars. 

1.4 The Problem of Complex Waveform of the Imput 
Voltage of DC and AC Machines 

This problem is of relatively recent orgin, since it is only in the 

last few years that controllers have been developed utilizing solid 

state components (thyristors or transistors) for supply, switching or 

chopping to yield a lower voltage level. However, some analytical 
(31) 

work was done as early as 1912 on the operation of DC series motors 

from an AC power supply. The performance of DC series motors when 

(32 33) 

controlled by pulses was extensively examined by Franklin ’ 
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with the aid of analytical tools. In this treatment the inductance and 
resistance of the machine was assumed to be constant, independent of 
frequency and saturation. A simple model of the machine was developed, 
which describes the operation under various speeds, pulse widths and 
periods. 

Such a model, although useful to describe the machine from a user's 
viewpoint does not appear to accurately predict the performance of a 
machine. Recently, DeWolf^ , and Hamilton et al^ } demonstrated that 

the inductance and apparrent resistance of a DC machine are functions of 
the current level and the frequency, and noted the need for a model 
which would encompass this phenomena. The external manifestation of the 
eddy currents is an increased resistance resulting from eddy current 
losses in both conductors and in iron portions in the magnetic circuits 
and the saturation effect is that of changing inductance as the currents 
change. 


1.5 The Continuous Modeling of the Electromagnetic 
Field - Incremental Inductances 

In both cases, of synchronous motor starting and of complex input 
waveforms, the problems associated with predicting, from design para¬ 
meters, the performance of the machine stem from the nonlinearities 
involved. These nonlinearities are caused by the eddy currents induced 
in the solid iron parts and the conductors and by saturation of the 
magnetic circuit. At any instant, the rate of change of induced volt¬ 
ages and of currents will depend on the input voltages and also on these 
nonlinearities to an extent that can be determined only by the knowledge 
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of the electromagnetic field at that instant. Incremental, or small 
signal, inductances can be calculated from the voltages induced in 
conductors due to small changes of currents. These inductances, together 
with terminal voltages, can be subsequently used to predict currents and 
voltages after a small time increment, and this process can be repeated 
and the machine performance calculated for any desired length of time. 

1.6 Statement of the Problem 

As has been discussed in the previous articles, the performance 
of a machine under continuously changing operating conditions poses a 
problem, both to the user and the designer which cannot be treated using 
traditional tools. The main causes of such an inability of conventional 
methods to deal with this problem are the presence of nonlinearities, 
namely the eddy currents in conductors and solid iron portions of the 
machine and the magnetic saturation of the iron of the magnetic circuit. 

In this dissertation this problem is addressed to as basically a 
transient problem from the viewpoint of the calculation of the electro¬ 
magnetic field. The method of time dependent finite elements is utilized 
to model the field in the time domain, and both mathematical models and 
computer techniques are developed to solve efficiently the resulting 
systems of equations and to calculate parameters which describe the 
performance of the machine, as currents in conductors, induced voltages 
and torque. 

While doing so, it becomes clear that electric machinery problems 
which appeared to be entirely unrelated, can be treated with the same 
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techniques. Two such transient problems are discussed in this disser¬ 
tation, the asynchronous starting of synchronous motors and the operat¬ 
ion of a chopper controlled DC machine, and essentially the same tech¬ 
niques are used on both to predict their performance. Such a uniform 
way of analyzing the operation of a machine gives a deeper understanding 
of the principles on which the theory of electrical machinery is based, 
and at the same time makes possible the calculation of such parameters 
as the proper shift of the brushes in a DC machine and the current dis¬ 
tribution in the damper bars of a synchronous motor. 
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2.0 THE MAGNETIC VECTOR POTENTIAL AND THE QUANTITIES 
OF THE ELECTROMAGNETIC FIELD 


The solution of the electromagnetic field inside an electrical 
machine is a difficult task, even when one is to use numerical tech¬ 
niques, because of the fact that the field is three dimensional and the 
geometry changes with time. In order to simplify the problem, the 
machine is assumed to be infinitely long and all the currents parallel 
to its axis. These assumptions mean that the field is two dimensional 
at any cross section of the machine, and simpler techniques can be used 
for its solution. The field in the end region cannot be considered two 
dimensional, at least at cross sections perpendicular to the axis of the 
machine. It has to be calculated separately, and the portion of the 
inductances due to it - the "end-turn inductances" - evaluated. The 
currents and voltages can be subsequently adjusted by taking these 
end-turn quantities into consideration. 

A similar problem is encountered when dealing with the eddy currents 
both in solid iron portions of the machine and in laminated steel. The 
eddy currents induced have to travel a distance longer than twice the 
length of the machine, since they have to form a loop by following a 
route transverse to the axis of the machine, as is shown in Figure 2-1. 
The conductivity assigned therefore to iron portions must not be the 
actual but a corrected one equal to: 


a 

corr 


2 1 + 


2 1 + 



a • 


(2-1) 
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2.1 The Magnetic Vector Potential 

From Maxwell's equations' - , in a domain without electric charges or 
polarized media, and neglecting displacement currents: 

V x B = p (J + V x M) (2-2) 

where p is the magnetic permeability of the medium, B is the magnetic 
flux density, J is true density, M the magnetic moment density and VxM 
the atomic magnetization currents. In the absence of permanent magnets 
this term can be dropped, and equation (2-2) becomes: 

V x B = p (J) (2-3) 

Also, from Maxwell equations: 

V • B = 0 (2-4) 

Equation (2-4) implies that B can be written as the curl of a function: 

B = V x A (2-5) 

where A is defined as the magnetic vector potential. In terms of A, 
equation (2-3) can be written as: 

V x (- V x A) = J 


(2-6) 
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The current density J may be either a known quantity, as is the case of 
applied currents in conductors, or unknown as is the case of eddy currents. 
Again, from Maxwell equations: 


V x E = - 


SB 

St 


From equations (2-7) and (2-5): 


(2-7) 


V x E = - 


(V x A) 


(2-8) 


Thus, E and - g- A differ by a function, the curl of which is zero. This 
function, therefore, can be written as a gradient of another function $: 


E 


A + grad $ 


(2-9) 


This last equation means that the actual electric field which can 
be measured by any suitable method, can be considered to be due to the 
rate of change of the magnetic field (as given by the first component of 
the right hand side of eq. 2-9), and to an applied electric scalar 
potential <t> (the second part of the right hand side of eq. 2-9). 

On the other hand, the current densities can be calculated from the 
strength of the electric field as: 

J = a E (2-10) 


where a is the conductivity of the medium. 
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From the above equations, it is clear that when instead of the true 
current densities, the applied electric field is known, equation (2-6) 
can take the form: 


V x (J V x A) = - a ( a | A - grad<t>) (2-11) 

9A 

The first term of the right hand side of eq. (2-11), -a^-, is the den¬ 
sity of the eddy currents, while the second part,-9 grad$, is the super¬ 
imposed current density due to a voltage, either applied, or present due 
to electrostatic phenomena. 

In the two dimensional case addressed here, equations (2-6) and 
(2-11) can be further simplified, since the only component of the current 
densities present is the z component, parallel to the axis of the machine 
This means that the z component of the flux density, B , vanishes from 
equation (2-4), and from (2-5): 


B 

x 


.2 A 
9y z 


(2-12) 


B 

y 


9x z 


(2-13) 


and, 


V x (- V x A) = - V(- V • A ) • k 

p H z 


(2-14) 


k being the unit vector in the z direction. 
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Equations 2-7 and 2-12 then become: 


V (- V A ) 

(J z' 


= J 


V (- V A ) 
|J z' 


= a 




84k 

8z J 


(2-15) 


(2-16) 


These two equations can be solved over the entire cross-section of the 
machine. Knowledge of the value of A^ makes it possible to calculate 
flux densities, induced voltages, eddy currents, losses, forces and 
torque. 


2.2 The Calculation of Induced Voltages 
and of Eddy Currents 

If a conductor is located in an electromagnetic field and is either 
moving with respect to the field, or the field quantities are changing 
with time, an electrical field results in the conductor, of magnitude: 

E = - ff + grad* (2-9) 

grad4> is the component of the field strength due to charges or external 

8A 

connections, whereas the induced voltage is - If the field solution 

ot, 

holds for a length £ the voltage along this length will be: 


V = £ (grad4>) z 


(2-17) 
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If that part of the machine which was called a conductor above, is 
a line in the iron of the rotor or stator parallel to the axis of the 
machine, this voltage must be zero. This can be shown with the aid of 
Figure 2-2. Assume A and A' to be lines symmetric with respect to the 
axis of a two pole machine, and parallel to it. Due to the symmetry the 
voltages, currents and magnetic vector potentials at A and A', and B and 
B' will be opposite. On the other hand, all points at every Cross 
section are connected together; therefore, they have the same potential. 
It follows that this voltage should be zero. This means that the eddy 
current density in solid portions of rotor stator are given by 


J eddy 



(2-18) 



Figure 2-2 The grad at symmetrical points. 

In an actual conductor, the value of grad4> is not identically zero. 
The field strength is given again by equation 2-9. When the current 
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density is constant over the cross section of the conductor, E is also 
constant. The grad4>, since it is externally applied, it is also uni¬ 
form over the cross section. 

In the case of solid conductors of relatively large cross sectional 
area, eddy currents can be induced, and the resulting skin effect can be 
pronounced. The skin effect can be treated in the same fashion as eddy 
currents. These eddy currents sum to zero and the total current remains 
the same as if there were no eddy currents, with an increased apparent 
resistance. 



Figure 2-3 The effect of eddy currents on the current 
distribution in a conductor 
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la order to account for these eddy currents, revisions are needed in the 
way in which the magnetic vector potential is calculated: 

Total current in a conductor, of cross-section S, is given by: 

I = fJ dS=-a/||dS+aJ grad $ dS (2-19) 

S S S 

I=-a a |/AdS+a grad $ S (2-20) 

S 

If the total current is known, the current density in a conductor 
takes the form: 

J = - a || + ° grad <l>*J=-°§£ + §- §- a f J A dS ( 2 “ 21 ) 
and the average induced voltage in the conductor is: 

V ind = - t a! I A dS ) • f < 2 - 22 > 

s 

These last equations can be simplified when numerical methods are 
used, as will be shown later. 

2.3 The Maxwell Stress Tensor 

In field theory it should be possible to calculate the net force on 
a given volume element within a magnetic field, using only the field 
conditions on the surface of the volume. This implies that the field is 
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a stress transmitting medium in the same sense that a string tying two 
weights together is the medium that transmits a force from one weight to 
the other. 

Following this procedure and considering that the force of a given 
volume is transmitted through the surface of this volume, the transmit¬ 
ting force can be formulated in terms of a quantity known as the Maxwell 

th 

stress tensor, T. The ap component of this sensor, T^p, is constituted 
th. 

so that the a component, dF^, of the force, dF, transmitted across a 
surface element, dS, and whose component is the P*"* 1 direction is dSp, is 
given by: 


F a = 


p> 


dS, 


(2-23) 


When the following restrictions are imposed for the medium through 
which the forces are transmitted: 

1. B is linear, i.e. permeability is not a function of the field. 

2. No permanent magnets are present. 

3. There is no magnetorestriction. 

then the stress tensor components take the form: 




(2-24) 


where 6ap is the Kronecker delta. 
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The matrix corresponding to this tensor is: 


T - r 


2 2 2 

\ (B - B - B ) 
x y z 


B B 
x y 


B B 
x y 


B B 
x z 


3s (B 2 - B 2 - B 2 ) B B 

y Z X y 2 


B B 
x z 


B B 

y z 


2 2 
% (B„ - B„ 


2 

B ) 
y 


(2-25) 


This tensor can be reduced to three components by transformation to 
principal axes. These axes are oriented so that one is parallel to the 
vector B and the two other are perpendicular to each other and the first 
axis: 



( 2 - 26 ) 


To illustrate, assume such a system of axes. One will be parallel 
to B and the other perpendicular to B, in the plane defining by B and 
the normal to the surface element as shown in Figure 2-4. It is seen 
from this figure that the magnetic field bisects the angle between the 
normal to the surface, and the direction of the resultant and force 
acting on the surface. 

In the application to electrical machines, the main concern is the 
torque acting on the rotor of the machine, rather than forces. A 
surface on which such a calculation can be carried out, is a cylinder 
lying in the air gap and enclosing the rotor as shown in Figure 2-4. 
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From the solution of the magnetic vector potential in the cross-section 
of the machine, the flux density can be calculated everywhere in the 
machine. In the two dimensional case, equation (2-23) and (2-25) become: 






(2-27) 


(2-28) 
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- « the rotor> M> can then be caicuuted by 

the surface of the cylinder by: 




M = 


d F 


(2-29) 



Figure 2-5 

V.S 


A surface on which the 
can be carried out. 


calculation of the torque 
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3.0 THE FINITE ELEMENT METHOD AS APPLIED TO THE 
ELECTROMAGNETIC FIELD CALCULATIONS 

Since 1970, considerable research has been devoted to applying the 

Finite Element Method to the solution of equations 2-15 and 2-16 at the 

cross-sections of electrical machines and to utilizing the result for 

the prediction of their performance. A general overview of the appli- 

(37) 

cations and of the theoretical background is given by Chari . What 
follows is a brief description of the general method, the mathematical 
analysis and the techniques that utilize the solution to calculate the 
field quantities and machine parameters. Both the time independent and 
time dependent solutions are discussed, since the electrical machinery 
problems are solved in the time domain. 

3.1 The Variational Method 

Consider a class of problems characterized by equations of the 

form: 

A(x) . u(x) = f(x) x € Q (3-1) 

where u(x) satisfies the boundary conditions: 


= 0 (or u=0) on the boundary 8Q 
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and A is a strongly elliptic 2 n< ^ order partial differential operator, 
i.e. it can be written as: • 


A(x) = E (-1)1“* D a a (x) D B 

|a|,|e|±2 


(3-2) 


where the operator D is defined as 


a 


n 


D a u(x) -= /d' r a„ 

3x^ 3 x 2 


3x 


n 


n 


(3-3) 


A bilinear form B(u,v) can be associated with A, where: 


B(u,v) - f E 




a Q (x) D a v D^u dx 
ctp 


. (3-4) 


It is possible now to replace the original problem of solving (3-1) 

« 

with the equivalent variational one of finding a u, such that: 


B(u,v) = (f ,u) ¥• u 6 H (fi) 

p 


(3-5) 




where the operation (\ ,.) denotes the (ft) inner product and H q (ft) 

is a Hilbert space of second order. 

The Galerkin method for the solution of (3-5) can be now described. 

Identify a subspace, (ft) of H o m (ft), spanned by a system of lin- 

G 

early independent functions {(^.(x). Each element of (ft) is of 
the form: 

0 

G . 

V(x) ^IC 1 <)>.00 
i=l 1 


(3-6) 



where C 1 i = 1, G are constants. 

The Galerkin approximation to u(x) is then a function 


U(x) = I A 1 <)>. (x) 


(3-7) 


which satisfies: 

B(U,V) = (f,v) v ve S h (Q) 


(3-8) 


Thus, by introducing (3.-7) in (3-8) a system of linear equations is 
obtained: 


2 K. . 

ij 



(3-9) 


where 





f. 
i 


(f 


(3-10) 


For the time dependent case, equation (3-1) becomes: 

9u(x,t) 

9t + A(x) • u(x,t) = f(x,t) (3-11) 

u (x,o) = u q (x) 


and the Galerkin approximation at every t becomes: 
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which satisfies: 


(§^—V) +B(U(t),V) = (f(t),V) 

te (o,T] 


(3-13) 


and the initial condition 


(U(o), V) = (u q ,v) 


(3-14) 


Again, by introducing (3-12) into (3-13) it is possible to obtain a 
system of first order differential equations: 


G 

I 

j=l 



A J (t) + K ± . A J (+)] f.(t) 


(3-15) 


where: 


K., = B (4>., (Jk) - • — - (3-16) 

f ± = (f, 0^) 


The system of linear equations (3-8) can be solved using various met¬ 
hods; an overview of these methods is given in the following chapter. 

S The system of first order differential equations, (3-15,) requires 

further consideration. 

In the case of linear coefficients, the Crank-Nicholson-Galerkin 
approximation can be used. To start, a partition P of the time interval 
(0,T) is introduced, composed of the set {t^, t^ 


, . .., t^} where 
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0 = t < < ... <t^ = T with t^ - t Q = At^ ; also the sequence 

{u n }^_Q is introduced to denote the value of the function U(x,t) at all 
the points of the partition P. Then, a family of finite-difference 
Galerkin approximations is introduced, associated with 6(0 < 0 < 1) 
which represents solutions to: 


( u n+1 -U D ,V) + (1-6) B (U n+1 ,V) +6 (U n ,V) = (f(+),V) (3-17) 

At 


Expansion of (3-18) using (3-14) yields: 


X[G.. + A-t(l-6)K. .]A n+1 = I (G. .-At-0-K. .)A n + At f. (3-18) 
I i] in j ij j i 


which for every n is a system of linear equations that can be solved 
with the same methods used for (3-9). In nonlinear cases, the differ¬ 
ential operator (3-2) becomes a function of the solution u. Denoting 


B(w;u,v) =j‘j a | f |g | l2 a a3 < x ’ w > D ° v dx 


(3-19) 


then, a two level predictor corrector version of the previously described 
Crank-Nicholson-Galerkin approximation is given by: 


" m+ L -— ' V + B( Vl <1+e) W m+1 + 7 - 0 


(3-20) 


Bmtl At Uj : • v + B( l (1+e) Vi + 7 (1 - 9) V 


I (1+9) Vl + 7 (1 - 0) "m- V > ’ 0 


(3-21) 


Such an approach requires the solution of two systems of linear alge¬ 
braic equations, for W m+ ^ and U at each time step. 


3.2 Computing Considerations and the Formulation 
of the System of Equations 


Throughout this investigation, linear triangular elements have been 
used, since they can be easily defined over any domain, and simple 
enough to be handled in a computer program. The details of the cal¬ 
culations associated with them are presented herein. 


3.2.1 The Local Matrix and Vector of a Linear Triangular Element 


Assume a triangle with vertices at (x^,y^), ^ 2 ^ 2 ), (x^y^). The 


linear interpolation functions 8^: 

((). = a. x + b. y + c. (3-19) 

T i i i J i 

should be 1 at the point i and zero at the other two points. From this 
consideration, the constants a., b., c. can be calculated as: 

’ l l ’ l 


a, b, c. 


x„ x~ 

111 


12 3 

a 2 b 2 c 2 

ii 

y i y 2 y 3 

a_ b„ c~ 


i i i 

3 3 3 


—i 


(3-22) 


From equation (2-5), the quantities B(<}>^,) are given by: 


B(<() 1 , (j).) =- / - 4 <J>. 4 <)>.=- — / a. a. d 

T 1 T j J p dx Y i dy p J 1 j s 



A . 
2 


(3-23) 


where: 
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Defining the integrals over an element: 

I 1 = Jx ds = (x 1 +x 2 +x 3 ) A/6 

i 2 = Jy ds = (yi + y 2 +y 3 } A/6 

13= Jxy ds = [(x^x,^) (y 1 +y 2 +y 3 ) + {x^+x^+x^ 
I^= Jx 2 ds = (x 2 +x 2 +x 2 +x 1 x 2 +x 2 x 3 +x 3 x 1 ) A/12 
i 5 = ; 2 ds = (v 2 1 + v 2 2 + y 2 3 + y l y 2 + y 2 Y 3 + Y 3 Y l ) A /12 

Then, 

(<b, , <j>.) = f(a,x+b.y+c . ) (a.x+b.y+e.) ds= 

Y J 1 1 1 J J J 

=[a.a.I, + b.b.I c + c.c.I A + (a.b. + a.b.) I 0 
1 j 4 l j 5 1 j 0 ij Ji 3 

+(c.a. + c.a.) I 1 + (b.c. + b.c.) I„] (3-26) 

i J J i 1 i J J i 2 

and for f constant over the element: 

f(a.I + b.I + c.) A/2 
i i i 


the local matrix has as its ij^ entry the term G„ + At(l-0)K^. 

local vector has as the i^ entry the term f. + (G. . - At0K. . )A? 

J x ij ij i 


,f) = /(+ b t y + c i )f 


If: 


G. . = - a (4>. ,()).) 
ij i J 

K.j = B (<J) i , 

= (f , «i) 


(3-25) 
)] A/24 


(3-27) 


(3-28) 


and the 
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3.2.2 Boundary Values and Negative Coupling 


It is a characteristic of electrical machines that at points having 
the same radius and differing in angle by one pole pitch, the flux den¬ 
sities are of equal value and opposite direction, and that the flux 
lines (equipotential lines) run parallel to the outer shell, which means 
that the whole shell surface is at a constant magnetic potential. By 
assigning zero values to the potential of the points of this surface, 
the potential at points shifted by one pole pitch assumes opposite 
values, with the result that calculation of the potential of points of 
the machine which lie outside an angle corresponding to a pole pitch is 
redundant. ' 



Figure 3-1. For the determination of the periodicity conditions. 
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In Figure 3-1 the value of A at 3Q^ and 8 Qg is zero. This results: 

A(r) 8 ft 3 = - A(r) 8 fi 2 = A(r) BQ^ (3-29) 

This case is treated by constructing compatible boundaries at 8 fl 2 
and 8 Q 3 , i.e. the grid is such that the nodes at the two boundaries 
correspond one-to-one and the corresponding ones have the same radius. 

Assume an element having distribution functions 4>|, 4> 2 > 4*3 and 
values at the nodes of Ap A^ and A^. The solution inside the element 
is then given by: 

A = Aj <t > 1 + A 2 <t > 2 + Ag (l)^ (3-30) 

The ijentry of the local matrix is: 

a^ = - o (4..,^) + (1-0) A t n B(<t> i , 4 k) (3-31) 

and the i *"* 1 entry of the local vector is „ . ----- 

S.-t i =(f > 4 ..) + [-a (4> i , 4>j)-9 At n B(4)., 4>j)l A n (3-32) 
When the value of the magnetic vector potential at one point, Ap 

I 

is the opposite of that at another point, A^ , equation (3-30) becomes: 

A = -a’4 » 1 +a 2 4 > 2 + a 3 4 > 3 = a^(-4) 1 ) + a 2 4 > 2 + a 3 4 > 3 


(3-33) 
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which means that a new shape function can b e used, and since: 


B (<!>., <b.) =-B (-4)., ,4» ) =-B (<{.., - 4 ^) . (3-34) 

($., 4>j)= -(-4> i , f.) =- (4> i » -tj) (3-35) 

the local matrix and vector become 




! f 

when A^ = -A^ and A ^ = -A^ , it follows easily from the previous 
considerations that: 



(3-37) 


(3-38) 
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3.2.3 Modeling the Airgap and Rotor Movement 

For the purposes of this study and the computer programs developed, 
the rotor and the stator are considered defined by coaxial cylindrical 
surfaces and include iron, copper, insulation and also air regions. 
Also, of course the airgap is defined by two coaxial cylindrical surfaces 
The air gap links two regions which move with respect to each other with 
the passage of time. For this reason, the whole domain - rotor, stator 
and air gap - cannot be discretized in a unique, time independent fashion 
If the air gap is excluded, the remaining two regions (rotor and stator), 
although moving in space, can have unique discretizations and global 
matrices, since they are independent of the position of the domain and 
its movement. 

In Figure 3-2 the discretization of the air gap is shown for three 
consecutive time instants. Assuming a discretization of the rotor and 
stator such that the number of points on their surfaces interfacing the 
gap are equal, SI is defined as the point on the stator slice modeled, 
which is the leftmost on the surface interfacing the gap, and R1 the 
corresponding point on the rotor surface. G1 is then defined as the 
midpoint between SI and Gl. Figure 3-2(a) then, is enough to define a 
discretization of the airgap for this initial position of rotor and 
stator. It must be noted that the points at the right hand side of the 
boundary of the slices, SN, GN and RN are to be coupled with SI, Gl and 


Rl. 
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(a) 


(b) 


(C) 


As the rotor moves, the elements of the air gap are reconstructed 
but are still defined by the same nodes. It is evident that the rec¬ 
tangles and elements of the air gap will be skewed for a period of time, 
but after the angle of point R2 has become closer to the angle of SI, 
the nodes defining the elements of the air gap must change, in the 
fashion shown in Figure 3-2(c). The negative coupling will remain as 
before, since the solution is sought for the initial slices of the rotor 
and stator. 

The initial relative position of the rotor and stator slice is 
shown in Figure 3-3(a). As in Figure 3-3(b), after a short time R2 
corresponds to SI and corresponding to -SI is point -R2. After some 
time lapse the slice of the rotor modeled is shifted completely outside 
the corresponding slice of the rotor. The correspondence of point will 
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be as shown in Figure 3-3(c) and it will be necessary to solve for the 
opposite of the values of the rotor surface. 

Such a movement results in changes in geometry of the gap, which in 
turn causes changes in the nodes and the local matrices of the gap 
elements. These changes do not markedly affect the global matrix; they 
only affect the entries related to the points at the surfaces of the 
rotor and stator, and the points in the air gap. 

Programming this movement in a computer is not as complicated a 
problem as it first appears. The treatment is based on a pointer, which 
shows which point on the rotor surface is corresponding to SI and 
whether it is coupled positively or negatively to the gap points. If at 
time zero, an angle of zero is assigned to all the points at the right 
side of the rotor and stator slices, and the angle of the points at the 
left hand side is a, and if 6 is the minimum angle between two consecu¬ 
tive points on the rotor surface, then the pointer moves down the rotor- 
air gap interface when the angle of point R1 exceeds a+6/2 (for clockwise 
movement) and changes sign every time it reaches the point ‘which was 
originally the rightmost. 

The local matrices of the airgap elements are formulated taking the 
negative couplings into account. Such a technique makes it unnecessary 
to solve in the whole rotor domain, which for a six pole machine would 
give a matrix prohibitively large to manipulate within a reasonable 


computer time. 
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3.2.4 The Formulation of the Global Matrix 

In every time step of the solution, the permeability of many of the 
materials and the conductivity of some of them may change. In order to 
avoid unnecessary calculations, the values of G^/a and M K„ in equations 
(3-28) must be stored, since they are independent of a and p. Each 3x3 
local matrix is symmetric, therefore, only six of its nine entries have 
to be stored. After every time step, p and possibly a are recalculated 
for each element, and subsequently the local matrices. The global 
matrix is then assembled, taking into consideration the negative coupl¬ 
ings of boundaries. The following example illustrates the way the 
matrix is assembled. 

Assume local matrices and vectors for the elements in Figure 3-4. 



Figure 3-4. A collection of finite elements; the node 
numbering outside the elements is global, inside is local 


\ 
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(3-39) 


The global matrix and vector become: 


a 33 

a 23 

0 

0 

0 

a 13 


a 32 

a 22 +b 33 

b 23 

0 

0 

a 12 +b 13 


0 

b 32 

b 22 +C 22 

C 12 

0 

b 12 +C 23 


0 

0 

C 21 

C ll +d 22 

d 12 

C 13 +d 23 


0 

0 

0 

d 21 

d ll 

d 13 


a 31 

a 21 +b 31 

b 21 +C 32 

C 31 +d 32 

d 13 

a ll +b ll +a 

33 +d 33 

(3-40) 
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V= 


a 2 + b 3 
b 2 + c 

C 1 + d 2 


a l + b l + C 3 + d 3 


(3-41) 


3.3 Numerical Calculation of Torque 

As discussed in the previous chapter, the torque can be computed in 
the air gap, half way between the rotor and the stator surfaces. For an 
element in that area with nodes at magnetic potential of A^, and A^, 
the potential over the surface can be given by: 

A= Aj 4*2 + A 2 ^2 + A 3 ^3 (3-42) 

(j>i = a i x + b i y + c i i = 1,2,3 (3-43) 


where a^, b^, c^ are given by eq. (3-22). From equations (3-42), (2-12) 
and (2-13), the flux density components inside the element can be cal¬ 
culated: 


A 1 b i + A 2 b 2 + A 3 b 3 


B y =- (A x a ! + A 2 a 2 + A 3 a 3^ 


B = 0 
z 


(3-44) 
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The stress tensor in Cartesian coordinates is: 



B - B 
x y 

2 B B 

y x 




(3-45) 


Then the forces applied on an elementary surface, dS, with components 

dS and dS are: 

x y 


dF '= T 

dS + T 

dS = 

1 

[(B - 

B ) 

dS + 2B B dS ]’ 

X 

xx x xy 

y 

2 % 

X 

y 

x x y y 

dF = T 

dS + T 

dS = 

1 

[2B B 

dS 

+ (B 2 - B 2 ) dS ] 

y 

yx x yy 

y 

2vi o 

y x 

X 

y x y 


(3-46) 



Figure 3-5 The calculation of torque in a machine 


The torque contributed by this elementary surface is: 


d M = r 


sin (}) • dF -r • cos 0 ♦ dF = ydF -xdF 
x y J x 


(3-47) 
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and the total torque is given by the integral over the whole rotor 
surface, thus: 


M = ( 


' ■ h~ # <B x - B y> ds 

c o •> 


+ B B dS 
x x y y 


+ 2^x B B dS + (fx (B 2 - B ) dS ] 
7 x y x 7 x y y 


(3-48) 


The rotor surface is represented by small linear segments, AS. The 

values of B and B are constant on each segment, therefore: 
x y 


£ [ (B 2 - B 2 ) / ydS + 2B B ( / ydS + /xdS ) 


o AS 6 c 


+ (B 2 - B 2 ) / xdS ] 

x y as y 


(3-49)> 


Equation (3-49) defines an algorithm by which the torque can be computed 
from the knowledge of the magnetic vector potential at the nodes of a 
grid, as calculated from the finite element) method: 

/ 

1. On an element, calculate B and B from (3-44) 

x y 

2. Calculate the torque contributed by the linear segment of 
the surface, AS, defined by the element 

3. Repeat steps 1 and 2 over the whole slice of the rotor that 
is modeled, adding the torque contributions. 

4. Multiply the result by the number of poles. 
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3.4 Calculation of the Permeability of Elements 

The permeability of each element is assumed to have a constant 
value over the element but it is also assumed to change (for magnetic 
materials) with a passage of time, depending on the flux density in the 
element. The simplest way to deal with this characteristic is to cal¬ 
culate the flux density from equations (3-44) and calculate the permea¬ 
bility from that and the saturation curve of the material. The solution 
(with the same force vector) is then repeated, and the permeabilities 
are recalculated. This procedure is repeated until the permeabilities 
calculated in the previous iteration are close enough to the new one. 
Various methods have been devised which accelerate this iterative tech¬ 
nique, the simplest of them being the underrelaxation method. In this, 
instead of the newly calculated value of permeability, M new > a linear 
combination is used, of M new and the previously calculated permeability 

M old‘ 


M = a p + (l-a)p ,, 
r r new 'old 

0 < a < 1 , usually 0.9 < a < 1 


(3-50) 


This method appears to be rather slow in the time independent problem, 
but gives good results for the time dependent problem and for iterative 
solution of the system of equations resulting from the finite element 
method, as Figures (3-6) and (3-7) show: 


No. of iterations 

in FEM solution Error in 
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Figure 3-6 Error in |J versus solution number 



Figure 3-7 


Number of iterations in conjugate gradient 
versus solution number 
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The case of laminated materials requires some special considerations 
The effective permeability of such a material will be less than the per¬ 
meability of iron, while the saturation level will he based on the flux 
density in the iron only, and not on the average flux density. The 
value of H will be the same in both iron and air. Assuming a stacking 
factor s, an iron permeability y, air permeability (J q , the average flux 
density B can be obtained using Figure 3-8: 


B = s • y * H + (1-s) ♦ y Q • H (3-51) 

B = [ s ♦ y + (1-s) y Q ] H 


From which the effective permeability will be 


H eff = t» • M ♦ (1-s) H 0 1 


(3-52) 


Conversely, when the average flux density is known from the solution of 
the electromagnetic field, as given in (3-42), the flux density in the 
iron, B. , can be calculated from (3-49) 

7 rnn 7 N ' 


B= s • B. + (1-s) y H 
iron o 

B= s B. + (1-s) y B. /y 
iron 'o iron r 

B. = B/ [ s + (1-s) y /y] 
iron r o r 


(3-53) 
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This value of the flux density in the iron, based on the previous value 

of the permeability, can be used to calculate a new value of the permea- 

/ 

bility 




Figure 3-8 For the calculation of the effective 
permeability of laminated materials. 

3.5 Calculation of Eddy Current Losses 

The losses due to eddy currents are of importance, particularly 
these in the iron portion of the machine. There, the current density 
is: 


CT 


3A 

at 


J= 


(2-18) 
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and it can be calculated for each element from equation (3-40) and• 
(3-41): 

3A 3A 3A 

J = - ff( *i3 r + ^2 ~3t + ^3 3t~" } 

= gr [<I> 1 (A^ +1 - A*) + 4> 2 (A n+1 - A n ) + <f> 3 (A- 1 - A")] (3-54) ■ 

n 


The losses in this element can be calculated from: 


dW= - J 2 dV dt 
cr 


(3-55) 


where dW is the Joule loss in an element of volume dV. For the case of 
two dimensional field, the losses in an element during the time step At^ 
are: 


AW = 


At 


SL f [A 


n 


n+1 .n . , . /.n+1 ,n . , , 

- A 1 ) h + (A 2 ' ~.. A 2 } ^2 + (A : 


n+1 


A 3^ 3 J Z dS 


(3 


(3-56)- 


where £ is the length of the machine and the integral extends over the 
cross section. 

In case of current carrying conductors, equation (2-18) becomes, 


J= 


<• 3 « 3*s 

0 ( 8t A z 3z^ 


(2-9) 


J= 


° si A z + 


I 

0 


(3-57) 


-56) 
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where I is the current flowing in the conductor due to external sources 
o 

Then equation (3-56) becomes: 


AW - [(Af 1 - Aj ) * 


1 + 


a 2 


+ (A^ +1 - A“ )<f> 3 ] + I Q } 2 dS 


(3-58) 
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4.0 SPARSE MATRIX ALGORITHMS APPLIED TO THE SOLUTION OF 
THE SYSTEM RESULTING FROM THE FINITE ELEMENT METHOD 

The system of linear equations derived from the application of the 
finite element method to an elliptical differential equation possesses 
certain properties, which can be utilized so that the solution is 
carried out in an efficient and accurate way. 

From the structure of the grid used in the finite element method, 
it is evident that every node is coupled directly only to the ones with 
which it is in the same element. In Figure 4-1, for example, node 1 is 
coupled <to nodes 2, 3, 4, 5 and 6, not to 7, 8, 9 or any other in the 
grid. As a result, the only possibly non-zero entries in the first row 
of the global matrix will be (1,1), (1,2), (1,3), (1,4), (1,5), and 
(1,6), no matter how large the total grid may be. Such a construction 
creates a global matrix of very large order (equal to the number of 
nodes), but also very sparse (very small percentage of non-zero entries).. 
Theoretical considerations also guarantee that this matrix is symmetric 
and positive definite. 

The computational techniques which can be used to solve this system 
of linear equations can be classified as either direct or iterative. 
Direct methods reach a solution after a fixed and predetermined number 
of computations, while in iterative methods, an initial guess solution 
is successively improved through a repetition of rather simple matrix- 
vector operations. These iterations are stopped when a specified 
accuracy is reached. 



Figure 4-1 A typical grid 

4.1 Storage Techniques 

The structure of the global matrix is of importance for the deter¬ 
mination of the method of storage. Two types of matrices can be encoun¬ 
tered; broadly classified as the banded and the generally sparse matrix 
(see Figure 4-2). In a banded matrix a pattern of non-zero terms occurs 
along the diagonal and both the storage scheme and the algorithm used 
for the solution can take advantage of this characteristic in order to 
decrease the storage area and increase the solution speed. 

The matrix is always stored in a compact fashion. This means that 
besides the non-zero entries, none, or only a few zero entries are stored 
Markers are used to point to the first and last non-zero elements in 




Banded matrix . Generally sparse matrix 

Figure 4-2 Types of matrices 

each row in the case of a variable bandwidth matrix, whereas for a 

A.’ 

constant bandwidth matrix, the number of vectors used is equal to half 
the bandwidth and each contains the entries of the matrix which lie on 
lines parallel to the diagonal. 

In the case of a generally sparse matrix, the storage scheme becomes 
more complicated. Three vectors are used to store the values of the 
matrix entries and the location of each within the matrix. One vector- 
is used to store the non-zero entries in the order that one would find 
them in the matrix while moving from left to right within a row, scann¬ 
ing rows from top to bottom. Another vector of integer values is used 
to point to the first entry of each row. A third integer vector indi¬ 
cates the column in which each entry lies. Since the matrix is sym¬ 
metric, only the upper triangular part must be stored. 

In Figure 4-3 an example is given of the storage scheme for a 
matrix 6x6. Vector A contains the values of the entries, vector JA 
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shows the row to which each entry belongs and vector IA points to the 
location of the first entry of each row. This storage method is very 
efficient from the storage space viewpoint, but the handling of the 
matrix requires many manipulations of the local index vectors, thus 
increasing the execution time. 



Figure 4-3 Storage of a generally sparse symmetric matrix 

4.2 Direct Methods of Solution 

A direct solution process can almost always be divided into three 
subprocesses: 

1. Triangularization of the matrix by factorization or elimination. 
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2. Back substitution on the right hand side. 

3. Iterative refinement to improve accuracy. 

Triangularization becomes more complicated, and the data handling 
more time consuming, as the form of the matrix departs from the constant 
bandwidth structure. It usually requires the storage on magnetic disk 
of the rows of the matrix - including some zero entries - and the con¬ 
secutive retrieval of a number of these rows at a time, so that the 
Gaussian elimination or the Choleski factorization can be carried out. 
The number of the zero entries that have to be stored also increases 
when the matrix is not banded. 

4.3 Iterative Methods of Solution 

These methods have the advantages over the direct methods in that 
their programming is simpler and that they efficiently utilize the 
sparsity of the matrix. The simplest iterative method, and the most 
widely used, is the Successive Over Relaxation and its variations. For 
a system of equations: 

A • u = b (4-1) 


with 


A= D-L-U 


(4-2) 
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where U vnd L are upper and lower triangular matrices respectively and D 
is a diagonal matrix, the following equation is applied: 

u n+1 = U n +iuD -1 ♦ (L*u n+1 + U u n - D u n + b) (4-3) 

where u n is the n*"* 1 approximation to the solution and uu is a real number 

between 1 and 2. This method is always convergent for symmetric positive 

definite matrices. The optimal selection of the parameter u) is of 

practical interest, since it should minimize the spectral radius of the 

(39) 

iteration matrix. For matrices possessing a special property A , the 
spectral radius is minimized for an u) given by 

U) = 2 / (1 + 7 l-p 2 ) (4-4) 


where |J is the spectral radius of the related Jacobi matrix D^(L+U). 
However, in practical problems, seldom is a matrix obtained having 
property A. An increased convergence rate can be achieved by Successive 
Line Over Relaxation (SLOR), where the computational molecule groups 
points along rows or columns of the mesh, thus forcing the system to 
have property A. When two rows or columns are included, the method is 
called the Successive Two Line Overrelaxation (S2L0R), and when both 
columns and rows are included in a form of a peripheral block, the 
method becomes the Successive Peripheral Overrelaxation (SPOR). These 
last three methods (SLOR, S2L0R and SPOR) require for their implementa¬ 
tion, a fairly regular grid, where the nodes are arranged in columns, 
rows, or circular patterns and the construction of the global matrix is 
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based on the shape of the grid. On the other hand, SOR, although rela¬ 
tively slow, does not require a regular grid and is often used for large 
irregular domains. 

(40 41) 

4.4 The P-condition number and the Preconditioning of a System * 

The asymptotic convergence of a system of linear equations with a 
symmetric positive definite coefficient matrix, has been shown to depend 
inversely on the P-condition number of the coefficient matrix, defined 
as the ratio of the maximum eigenvalue to the minimum eigenvalue of a 
positive definite matrix. Thus, to improve convergence, it is required 
that the P-condition number of the matrix be very small. When this does 
not happen it is desirable to "prepare" the system so that its condition 
number is minimized. This means that instead of solving the system of 
equations (4-1), it would be easier to solve system: 

Q • A • u = Q • b (4-5) 

where Q is a non-singular matrix and the condition number of the result¬ 
ing system is less than that of the original system. One such matrix is 
(I-uuL)^, and the system (4-5) can be rewritten as: 


[(I - uL) 1 • A • (I - ojL T ) 1 ] [(I - o)L T ) • u] = (I - u>L) 1 b 


(4-6) 
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Denoting: 


T 

v= (I-ujL ) *u 

d= (I-ioL)' 1 -b (4-7) 

B= (I-toL)" 1 *A*(I-uiL T )“ 1 

then B can be shown to be also a symmetric positive definite matrix. A 
minimum of the value of the P-condition number exists for u) between 1 
and 2, whereas for u) = 0, the system reverts back to its original form. 

• (42 43^ 

4.5 The Conjugate Gradient Method 1 ’ ' 

In gradient methods in general, instead of solving the system of 
equations (4-1) as is, a variational method is derived based on the 
observation that the solution of (4-1) is the same vector that maximizes 
the quadratic functional: 

F(u)= u^ • A ♦ u + u^ • b (4-8) 

This functional defines a family of similar ellipsoids in a Euclidean 
n-space with center at A V The solution, therefore, can be attempted 
by proceeding to this center, using a sequence of displacements, achiev¬ 
ing after each, a better approximation of the maximum: 

(n+l)_ (n) (n) 

u — u +fft) 
n r 


(4-9) 
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where is the approximation, p^ is the direction followed in 

the n^ step and is a real number selected so that F(u^ n+ ^) is 
minimized. At the n^ step a residual, r^ n \ can be calculated as: 

r (r ^=b-A • u (n) (4-10) 

The way that the direction vector p^ n ^ is determined defines the 
method used. The most obvious choice is the direction normal to the 
surface of the ellipsoid, i.e. p^ n ^ is chosen to be -r^ n ^. This tech¬ 
nique defines the steepest descent method and it is illustrated for a 
two dimensional case in Figure 4-4. 



Figure 4-4 Steepest descent method. 

This is an iterative technique, and the rate of convergence and 

\ 

number of iterations depend on the form of the equations and the accur¬ 
acy desired. 

A better method results when advantage is taken of the property of 
quadratic functionals, that the center of the ellipsoid lies on a plane 
conjugate to a given cord. 
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4.5.1 Conjugate Directions 

Given a symmetric matrix A, two vectors, and d£» are A-orthogonal 
or conjugate with respect to A if: 



(4-11) 


If a set of vectors d , ...,d, are A orthogonal, then these vectors 

O K 

JL 

A 

are linearly independent. A solution then, u , of Au=b, or equivalently 

T T 

that maximizes (-1/2 u Au-u b), can be expanded in terms of n vectors d Q 
... d^ i (A is a nxn positive definite matrix) as: 


* 

u = ct d +... 
o o 



d 

n- 


1 


(4-12) 


Multiplying by A and taking the scalar product with d^, yields as: 



The last equation indicates that the coefficients CT of equation (4-12) 

JL 

r\ 

can be evaluated without knowing u , which, however can then be calcu¬ 
lated from (4-12). Thus: 


4 n-l d! ■ b 

u = Z - d ± (4-14) 

i=0 d. • A • d. 

i l 
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The basic idea imbedded in (4-14) is that when a set of d 's is 

l 

selected, the unknown vector u can be calculated from it. This 

* 

expansion of u can be considered as the result of a process of n steps, 

where at the i^ step, a.d. is added. 

ii 

Viewing the procedure in this way, the conjugate direction method 
can be formulated in n steps; the k^ 1 of which is described by equation 
(4-15): 


u 


k+1 


= u k + 


°k d k 


(4-15) 


where: 



( 4 - 16 ) 


and 


r R = A • u fe - b .■— — - ~ " • - - (4-17)— 


This holds true because the difference between an initial guess, 
u , and the solution, u , can be written as: 


* 

u - u = a d + 
o o o 


a d 
n-1 n-1 


(4-18) 
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Multiplying by A and taking the scalar product with d^, as in (4-13) 
yields: 


T -k 

d k • A (u - u q ) 

°k ,T . . 

d k * A d k 


(4-19) 


Following the iterative process (4-15) from u up to u, , one obtains: 

O K 


U, - U = ct d + a. d, + ... + o, , d. . 
k o o o 11 k-1 k-1 


(4-20) 


and because of the A-orthogonality of the vectors d^: 


d k A • (u k - u o> = 0 


(4-21) 


Substituting (4-21) in (4-19) produces: 


T * T 

d ■ A • (u - u ) r k d k 



(4-22) 


An important property of the conjugate direction method is that in each 
step the functional (4-8) is maximized on the line: 


u 


= Vi + 


ad 


k+1 ’ 


(4-23) 
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where a is a real number, as well as in the linear space u q + B^, where 
is a linear combination of d Q , ..., d^_j meaning that the remainder, 
r, is perpendicular to the space B, as shown in Figure 4-5- 



—-Figure 4-5 Conjugate direction method ■■ 

4.5.2 The Conjugate Gradient Algorithm 

The conjugate gradient method is obtained from the conjugate 
direction method by selecting the successive direction vectors as a 
conjugate version of successive gradien£s which are obtained as the 
method progresses. Thus, the directions are not specified before the 
solution starts; they are determined sequentially at each step of the 
iteration. At step k one evaluates the current negative gradient vector 
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The last step of the algorithm is of particular importance when the 
coefficient matrix is of large order. From the theory of the conjugate 
direction method it is clear that the solution will be reached after a 
number of iterations equal to the order of the matrix. But, because the 
method is based on gradients, the algorithm makes good uniform progress 
towards the solution at every step, in contrast to the case of arbitrary 
conjugate directions, where the progress may be slight until the final 
few steps. In the conjugate gradient method, therefore, the progress is 
quite good at the first steps, whereas the last steps improve the solu¬ 
tion only slightly, making the last iterations unnecessary when the 
desirable accuracy has already been achieved. A measure of the distance 
of an approximation from the exact solution is the gradient. When the 
solution is obtained, the gradient vanishes and a norm of it can be used 
to evaluate its length and be compared to the desired accuracy. 

A very important advantage of the conjugate gradient method is the 
very simple formulae which are used to determine the new direction 
vector. This simplicity makes the method only slightly more complicated 
than the steepest descent method. 

4.6 The Preconditioned Version of the 
Conjugate Gradient Method^^’^^ 

When the matrix A in the functional (4-8), is ill conditioned, the 
ellipsoids, as in Figure (4-4), become long and thin, and thus the 
normals to the surface become coplanar. Preconditioning reshapes these 
ellipsoids, so that they take a less prohibitive form, thus enabling the 



■' ,N3 j (-*. 
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solution algorithm to lead more rapidly to a good approximation. An 
example is given here of the effect of preconditioning on the form of 
the quadratic functionals and the shape of the corresponding ellipsoids: 
Assume the system of equations: 


- 

- 


* — 



1 

-0.98 


X 1 


1 

-0.98 

1 


X 2 


1 








(4-26) . 


The solution of this system is equivalent to maximizing the 
functional: 


t x 



x 2 ] 


■ "l " -0.98’ 


' X l‘ 

[ Xf x 2 ] 

l‘ 

-0.98 1 


x 2 


1 

__ -- - 




- 


(4-27) 


2 , 2 

X 1 ~ 1-96 x 1 x 2 + x 2 - 2 X;L - 2 x 2 = 2c (4-28) 

The ellipsoids corresponding to equation (4-28) are shown in Figure 
4.6(a). 

When preconditioning is applied, the system of eq. (4-26) becomes: 


1 0 


1 -0.98 


1 .98u 


1 0 


1 

. 98to 1. 


-0.98 1 


0 1 


. 98oj 1 


1 


(4-29) 
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and the corresponding functional is: 


O 


x^+2 0.98(l-u)x x 2 +(0.98 2 -2*0.98a)+l)x^-2x 1 -2(1-0.98a))x 2 =2c 


(4-30) 


which has as corresponding ellipsoids for U) = 1.05 the ones in Figure 
4-6(b), and.for u) = 1 the ellipsoids in Figure 4-6(c) which are circles. 



Figure 4-6 A family of ellipsoids before and after preconditioning 

When preconditioning is used, the arithmetic operations per iterat¬ 
ion are increased drastically, but, due to the faster convergence, the 
overall speed is improved. The ^matrix B = (J-uuL) A(I-U)L ) is never 
calculated-, since such a transformation would destroy the sparseness of 
the matrix and make storage requirements prohibitively large. Instead, 
back-and-forward substitution is used every time the matrix B is to be 
multiplied by a vector. The algorithm described in section 4.3.2 re¬ 
mains essentially the same, while certain manipulations are added to 
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account for the preconditioning. The system has to be prepared so that 
the entries along the diagonal become one. The final algorithm takes 
the following form: 


1. Prepare the system by multiplying both sides with matrix G: 

m 

G.. = 0 for i ^ j 
ij 



Using for simplicity the same symbols for system (4-27) as 
for (4-1), but noting that the diagonal is now the unity 
matrix: 

2. Calculate v = (I-uiL^) u 

o o 

3. Calculate d = (I-U)L)' ^ b 

4. Calculate the remainder of the system: 

B • v = d as: 

v = d - B . v = d - (I-uiL) 1 . A . (I-u>L T ) V 
o o o 

5. Find the first conjugate direction as: 
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jO o 
d = -r 

and store it in the place of d. 

6. Calculate a = r • d/(d^ • B • d) 

as a = r • d/[d • (I-cL) 1 • A * (I-tuL 1 )" 1 • d] 

7. Calculate the new approximation to the solution: 
v: = v + a ♦ d 

8. Calculate the new residual 

r: = r -a • B • d = to ♦ (I-U)L)^ • A • (I-u)L^) ^ * d 

It must be noted here that B • d has already been calculated 

in step 6. 

9. T = -r * d/d^ ♦ B • d 

10. The new direction vector, d, Can be calculated as 

d = r + X • d 

11. The distance of the current approximation to the solution, V, 
from the exact one can be checked by finding the norm of the 
remainder. If it is larger than the desired accuracy, the 
algorithm is repeated from step one. Otherwise, the solution, 
u, is calculated in the next step. 

12. u = (1 - ujL^)^ • v and u: = G • u 

4.6.1 Computer Considerations 

The preconditioned conjugate gradient method requires the storing 
of four vectors, v, r, d and B • d. Involved in the preparation of the 
matrix are two matrix-matrix multiplications and two scalar-vector 
products for the determination of A, u and B in step 1. In step 2, a 
matrix-vector product is involved, and in step 3 one back substitution. 
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At the exit (step 12), a back substitution is required, as well as a 
scalar-vector product. In each iteration, a back substitution, two 
inner products and three scalar-vector products are involved in compu¬ 
tational effort. 

As later discussed, often solutions are required using the same 
stiffness matrix A and entirely different force vectors b. In these 
cases the initial multiplication is not repeated; instead the values of 
G are stored in the diagonal of the matrix A, while the adjusted values 
of L are stored in the lower triangle of the same matrix. Therefore, 
only vectors d and u have to be multiplied by G and G ^ respectively. 

When the saturation level of the iron parts in the machine is to be 
determined, the linear system is solved repeatedly with the entries of 
the matrix changing slightly, and the force vector remaining constant. 
These conditions cause the'solution vector of the previous repetition an 
excellent guess for the following. This policy drastically decreases 
the number of iterations required for each solution. 

4.6.2 Determination of the Solution Technique for the Continuous 
Modeling of the Electromagnetic Field. 

The fact that the solution at each time step is close to the solu¬ 
tion of the previous time step makes iterative solutions more attractive 
than direct ones, since the number of iterations is decreased by the 
introduction of a good guess solution. The storage requirements are 
less restrictive in the iterative techniques, since no space needs to be 
allowed for fillings during the triangularization. This is of special 
importance in this case, since the grid is quite irregular and the 
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number of fillings would be large. The simplicity of programming is 
also a significant factor in deciding in favor of an iterative scheme. 

Among the iterative techniques which are commonly used, SOR, is the 
simplest and usually the first choice. The matrix from the finite 
element method, though, does not possess the property A, and it is 
difficult to find the optimal parameter U). When the grid is regular, 
the matrix takes a very convenient form, making the variations of SOR 
very attractive. However, the limitation of a regular grid deprives the 
finite element method of the ability to utilize a fine mesh only where 
it is necessary,'while using a coarser grid where the solution does not 
vary greatly in space. 

The electromagnetic problem has the peculiarity that the entries of 
the local matrices of neighboring elements can differ by two or three 
orders of magnitude when the material is magnetic in one element and 
non-magnetic in the other since these entries are inversely proportional 
to the magnetic permeability. Such vast differences in numbers which 
are added together in the global matrix, make the global matrix ill 
conditioned. This characteristic makes it necessary, when SOR is used, 
to alternate between SLOR and SPOR^ . Preconditioning appears to be 

an elegant way of rectifying or improving this situation. 

(43) 

Reid has shown that the conjugate gradient method in its origi¬ 
nal form (uj=0) is a powerful algorithm, and with the further application 
of preconditioning, both in this work and other it has been proven 

both fast and convenient, especially in the case of time dependent 
solution of electromagnetic fields. 
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5.0 THE OPERATION OF CHOPPER CONTROLLED DC SERIES MOTORS 

Both experimental work and problems noted with the application of 
DC motors have made it clear that there exists a discrepancy between the 
analytically calculated inductance of DC machines and that observed when 
the motors are driven by a chopper. Also, recent experimental work 
(34,35) k ag s }j 0wn that a single value of inductance cannot be assigned 
to the machine; rather, both apparent resistance and inductance vary 
greatly with the frequency of the current which is used to measure them 
and with the value of the DC current on which the AC current is super¬ 
imposed. These effects are attributed to the combined phenomena of 
magnetic saturation of the iron parts of the machine, eddy currents in 
both solid and laminated iron and the eddy currents and skin effect in 
solid conductors in the machine. 

These phenomena cannot be adequately described by a small set of 
equations except in an incomplete and empirical form. Instead, the 
solution of the field inside the machine should be sought at every 
instant, and the values of the inductances and other parameters should 
be calculated from it and be used to calculate the quantities of inte¬ 
rest to the designer and analyst. 

When the inductances are known as a function of input and time, the 
currents can be calculated as well. The theory of the two axes cannot 
be applied for the calculation of the torque since voltages are induced 
in parts of rotor and stator which are not part of the windings. Instead, 
the theory of the Maxwell stress tensor is used, as applied to a dis¬ 
cretized grid (section 3.3). The radial flux density at the air gap can 
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be calculated in order to determine the optimal location of the brushes 
for linear commutation and the electromagnetic field can be solved based 
on that position of the brushes. 

5.1 Calculation of Inductances and an Algorithm for the Calculation of 

Currents and Torque 

In the absence of eddy currents, the small signal, or incremental 
inductance around an operating point of the machine, can be found by 
first solving the electromagnetic field for the actual current at that 
operating point, obtaining in this way both the values of the magnetic 
vector potential inside the machine and the permeabilities of all the 
elements which correspond to that current level. Then, the current 
changed by a small increment, the field calculated again and the flux 
linkage (with all windings) calculated. The ratio of these flux link¬ 
ages to the change in current gives the incremental self and mutual 
inductances, on which the calculation of the actual change of current 
can be based, for a small time step, provided the changes in time and 
current are relatively small (so that the slope of the saturation curve 
remains constant for this change), the change of currents calculated on 
the basis of this incremental inductance will be accurate. Thus 


" L*. 


AA 


1 “incr AI. 

incr 


(5-i>; 


AI = V/ 


mcr 

At 


(5-2) 


where A^ ncr is a small change of current, AA the corresponding change 
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in flux linkages, V the induced voltage, and AI the actual change in 
current. 

However, when eddy currents are present and of a magnitude which 
greatly affects the performance of the machine, the current change, AI, 
cannot be based on the same formulae described above. 

The reason is that in this case, the machine, in addition to the 
current carrying winding, has an infinite number of other windings 
representing the paths that the eddy currents can follow. These 
windings are not open circuited; rather each is short circuited through 
the resistance corresponding to that current path. Each of these 
windings has a self inductance and a mutual inductance with every other 
current path, as well as with the real winding of the machine. Even in 
the case when all these current paths are replaced by a finite number of 
circuits, in order to make the problem less formidable, the number of 
independent windings required is still too large to give an inductance 
matrix which can be handled efficiently. 

If, in some way, the current change had been known (either measured 
or calculated), then a value of inductance could be obtained, depending 
not only on the saturation level, but also on the actual change of 
current, AI, and the time step, At: 


L. 

incr 


_V 

AI 

At 


(5-3) 


Equation (5-3) can be used in an iterative way in order to calcu¬ 
late the change of current from the finite element solution. Initially 
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a "guess" change of current is applied and the field is solved, assuming 
this change to be the exact one for the time step examined. From the 
solution, the voltage induced in the winding is calculated and a value 
of the incremental inductance is obtained from (6-3). It is then sub¬ 
stituted in equation (5-2) to yield a corrected value of the change in 
current. The process is then repeated to yield a closer approximation 
to the correct incremental inductance. The iterations are stopped when 
the values of inductance obtained from two consecutive iterations are 
sufficiently close. 

This process can be very slow, with the values of inductances 
approaching asymptotically the final value, but for the calculations of 
each, the field has to be solved. When the field is solved in an itera¬ 
tive way, however, the previous solution is a good approximation to the 
next, decreasing the amount of the iterations needed for convergence. 

The number of iterations for the inductance calculations can be decreased 

(44) 

if the Aitken extrapolation is used to accelerate the convergence. 

5.1.1 The Aitken Extrapolation 

When the differences between successive values of the incremental 
inductance form a geometric progression, of ratio Of: 



a 


(5—4) 


where L., L„ L^, L, are the values of the incremental inductance calcu- 
12 3 4 

lated in the respective iteration; then the final value of L will be: 



Application of this technique decreases the number of the itera¬ 
tions for the inductance calculation. There are three solutions of the 
field required for the calculation of a (two of which require fewer 
iterations in the preconditioned gradient method), and one solution 
after the final value of the incremental inductance is calculated, so 
that the accuracy of it is checked and the field is solved for the final 
value of current so that eddy losses and torque can be calculated. 

5.1.2 The Rotational Voltage 

When the flux density is known in the machine, then the rotational 
voltage can be calculated for each conductor as: 



where is the radial component of the flux, r the distance of the con¬ 
ductor from the axis of the machine, £ the length of the conductor and a 
the angular velocity of the conductor. When the magnetic vector poten¬ 
tial is known, B can be obtained from it as: 
r 
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and the voltage induced in all the conductors will be 



3 




n 

L 

i=l 




3A 

0) • i 

3<f> 


i 


(5-8) 


where n is the number of conductors in the armature and 5 - 5 -[the deriva- 

9<t> \l 

tive of the magnetic vector potential at the conductor. If the distri-. 
bution of the winding is assumed linear over the armature,then: 


V = - n a) 


f n M 

3<t> 


d<j> = - n 


0 ) 


(A„ - A,) 


(5-9) 


where 1 and n are the first and last conductors of the uniform distrib¬ 
ution of the winding and c()^ and <{> n their respective angles. 

5.1.3 The Voltage Input and External Connections of the Machine 

The previously described method of handling the operation of a DC 
machine can be applied to any case where the waveform of the input •••• 
voltage is complex. Here the case of a series traction motor is 
examined when driven from a set of batteries through a chopper 
controller of relatively low frequency (100-400 Hz). The connections of 
the machine are shown for this case in Figure 5-1. 

The armature terminals are shunted by a reversely biased 
(freewheeling) diode, D2, and the armature plus series field shunted by 
D1. The first freewheeling diode, D2, does not carry any current in the 

• 1 


i 




case that is being examined here, whereas the second diode, D-^, conducts 
during the time that the chopper is not conducting. 


Chopper 



Figure 5-1 The external connections to a chopper controlled DC series 
motor. 

The open circuited voltage of the batteries is known and is assumed 
to remain constant during the simulation, i.e. the batteries are not 
supposed to discharge significantly. The battery impedance cannot be 
generally considered a constant, but for the purposes of this simulation 
it was ignored, since its value does not greatly affect the performance 
of the machine. 

5.1.4 The Algorithm 

The algorithm that was used for the calculation of the currents in 
a DC machine makes repeated use of the procedures described previously^. 
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in this section. After every time step the value of the resulting 
current is injected into both the field and the armature windings and 
the field is calculated. From this solution the Maxwell stress tensor 
can be calculated in the air gap, and the torque computed. Then the 
procedure of calculating the inductance and change of current is repeated 
The iteration stops after as many cycles as are needed for the current 
waveform to come to a quasi-steady state. The algorithm used is: 

1. Start with an assumed current increment, I. , and an old 

incr’ 


value of current I 


old' 


I: = I 


old 


(5-10) 


2. Solve the electromagnetic field based on current in the arma¬ 
ture and field, equal to I. 

3. Calculate the voltage, V.^, induced in the conductors due to 
the change AI in the time interval At and compute the 
incremental inductance: 


L. 

incr 


V / (—) 
ind' L At ; 


(5-11) 


4. Calculate a new change in current based on the new incremental 
inductance and a rotational voltage constant, K, calculated in 
the previous time step 
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where L . is the external and end turn inductance and R the 
ext 

resistivity of the winding 

5. Repeat from steps 3 three times for I - 1^^ + AI; calculate 
three values of incremental inductance and the final value 
from Equations (5-9) and (5-5) (Aitken Extrapolation) 

6 . Repeat steps 2, 3 and 4 for the new value of 

7. Calculate the air gap torque and the eddy current losses in 
the iron. 

8 . Calculate from (5-9) the rotational voltage, V rQt and the 
constant K as: 

K = V Jl 
rot 

9. Repeat from step one using, as the new value of I and 

as the initial guess of AI, that one calculated last. 

5.2 Application to an Electric Vehicle Motor 

The algorithm of Section 5.1.3 was applied to a motor which was de¬ 
signed for operation in electric vehicles, operated from an 84 volt 
battery supply, utilizing a chopper with maximum frequency of operation 
of 400 Hz and a minimum time on of 1 ms to control the average voltage 
applied to the motor. The machine was connected to a dynamometer instru¬ 
mented to measure the average current, output torque, and the revolutions 
per minute and loaded. Oscilloscope pictures were taken of the current 
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and voltage waveforms for various values of load, speed, chopper period 
and chopper conducting time. 

5.2.1 Characteristics of the Motor 

Prior to testing, the motor was dismantled and its dimensions 
measured and recorded. After reassembly, the friction and windage 
losses were measured and an interpolation technique applied to the 
results yielding the curve shown in Figure 5-2. Table 5-1 lists the 
main parameters of the motor. 

The armature was wave wound and the pole windings series connected. 
The brush holder was moveable, allowing the brushes to be offset up to 
30 electrical degrees from their neutral axis in either direction. The 
computer simulations presented here were with the two field windings 
connected in series and the brushes at the neutral position. 

Figures 5-3 to 5-5 show the significant dimensions of the rotor, 
stator, field winding, armature slot and armature conductors. 
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Figure 5-2 Friction and windage losses 
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Table 5-1 Parameters of the motor 



Figure 5-3 Dimensions of rotor and stator (in millimeters) 
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5.2.2 Approximations and Assumptions 

The total number of slots in the armature of the machine is 27, 
making the number of slots per pole equal to 6.75. Since this number is 
not an integer, an accurate representation of the machine would be 
obtained only if the complete cross section was modeled. The accuracy 
would not be significantly decreased if only one pole was modeled using 
an integer number of slots, and keeping the ratio of the length of slot 
to tooth the same as in the actual machine. 

In doing so, the number of Ampere-turns of the armature would be 
increased if the actual current of the machine was to be injected in the 
slot conductors. In order to simulate the field correctly and to main¬ 
tain in the model the same magnetomotive force as in the real machine, 
the current in the armature winding should be multiplied by the ratio of 
the actual number of slots per pole to the number of slots in the model; 
in this case by 6.75/7.00. This would ensure that the calculated values 
of the magnetic vector potential values and flux densities would be 
correct. 

After the field is calculated, the voltages induced in the armature 
due to the field winding (and its self inductance) can be computed. 
These voltages will be the sum of the voltages induced in the conductors. 
A correction is needed here also, since the number of conductors in the 
model is larger than in the actual case. The corrected voltage induced 
in the armature, due either to rotation or self inductance should be the 
voltage calculated from the model, multiplied by 6.75/7.00. 
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Commutation is a rather complicated phenomenon and it should be 
treated as such when it is studied alone. In that case the self induct¬ 
ances of the coils undergoing commutation are to be calculated, as well 
as their mutual inductance with each other and with the armature. Then, 
a model (consisting of differential equations) can be constructed based 
on these inductances and the resistivities of the various current paths. 
The problem is even more complicated by the fact that the brush resis¬ 
tivity is both nonlinear and dependent upon the brush pressure. On the 
other hand, if high resistivity commutation is assumed, its effect on 
the overall performance of the machine becomes minimal and can be 
neglected; i.e. the commutation assumed linear in the calculation of the 
currents of the coils undergoing commutation. 

The conductor portions of the machine require more considerations. 
As described in section 3.2, when no eddy currents are anticipated, the 
real current densities can be assigned to the elements in the conductors 
and their conductivity assumed zero for the eddy currents. This is the 
case when a conductor is composed of many fine strands. When eddy 
currents and .skin effects are anticipated, as is the case of solid ... 
conductors of relatively large cross sections, there are two possible 
ways to take their effect into account. One is to calculate the eddy 
currents in each individual conductor, by defining a relatively fine 
grid and applying equation (2-21) to all such conductors. The alterna¬ 
tive procedure is to assume that all slots are going to behave similarly 
and solve, in a separate grid, the problem for only one slot at the same 
time that the main grid is solved, and from that solution calculate the 
changes of apparent resistance. In this investigation, the skin effect 
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was neglected and the apparent resistance of the windings used was the 
resistance that was measured at 100 Hz. 

Figure 5-6 shows the discretization of one pole pitch. As can be 
noted, the discretization of the air gap is not shown. It is rather 
generated by the program as described in section 3.2.3. All conductors 
in the armature carry current in the same direction. Shifting of the 
brushes is accomplished by rotating the armature grid accordingly and 
letting the program define the airgap grid. 

5.2.3 Results of the Computer Simulation 

The program was run for various cases for which experimental results 
had been obtained previously. In Figure 5-7 the flux lines and the 
outline of the material boundaries are shown for the case of solid frame 
stator, where eddy currents are allowed to flow. In Figure 5-8 the same 
cross section is solved with the exception that the frame is assumed 
laminated. This is not the case of the motor discussed here, but the 
flux density distribution was examined because of interest in the 
characteristics between motors with laminated and those with solid 
frames. 

In Figures 5-9, 5-10 and 5-11 the results of the simulation are 
presented for a sample case, where only ten points were used for one 
cycle. These results are presented here for comparison with experi¬ 


mental results. 




Figure 5-6 The grid for the electric vehicle motor 
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Figure 5-9 Voltage waveform 

In Table 5-2 both observed test results and results of the simu¬ 
lation are given for a number of operating modes of the machine. The 
coarseness of the grid used does not seem to greatly affect the induct¬ 
ance and currents, but it obviously makes the eddy current calculation 
—inaccurate, as can be observed from Figure 5-7. It also affects the 
accuracy of the torque as calculated from the Maxwell stress tensor. 

The discrepancies between actual test data and results from the 
computer simulation can be attributed to the lack of knowledge of cer¬ 
tain material parameters. The torque and current calculations depend 
not only on the coarseness of the grid, but they are also very sensitive 
to the resistivity and permeability of the frame. Typical values were 
used in the simulation from which the results of Table 5-2 were obtained, 
although significantly better results were obtained when other values of 




CURRENT i TORQUE 

0 80.00 120.00; I. J3.00 20.00 40.00 60.00 







<D 






> 






cd 


LO 




H 

rH 

o 

Os 

O 

CO 


Os 

CM 

vO 


CM 

H 

• 

• 

• 

• 

• 

< 

o 

tH 

o 

o 

rH 



CM 

vO 

CO 

VO 

<u 

• 

• 

• 

• 

• 

> <1 

CM 

LO 

LO 


00 

cd 

CM 

i—I 

os 

CM 

vO 

H 

CM 

1—1 

rH 

rH 


/*N 

vO 

CO 

Os 

CM 

CO 

p- e 

vO 

<r 

vO 

o 

*sf 

cd 

• 

• 

• 

• 

• 

m a 


t—H 

Os 

CM 

CM 

s w 

00 

CO 

VO 

CO 

rH 



Cd cd( rH 
bO U| 

[J *K ' 



VO 

00 

CO 

i'g 

o 

CO 


vO 


t — 1 

• 

• 

• 


• 

o 

rH 

o 

o 

rH 



CM 6 


cm m 

o cm 

rH r^. 


O O 

CM CM 

rH oo 

pH t—I 



Table 5.2 Test and simulation results. 






















90 


resistivity and permeability of the frame were used. Also the effect of 
the batteries was completely ignored, although battery characteristics 
affect markedly the performance of the machine. 

5.2.4 Airgap Radial Flux - Brush Positioning 

The magnetic field due to the field winding in a DC machine is 
symmetric around a pole center line and almost uniform under the pole 
face because of the constant reluctance of the air gap (provided that 
the slot effect is neglected). Outside the pole face, the radial flux 
density drops rapidly, becoming zero in the middle of the interpolar 
space. When the brushes are at the geometrical neutral axis, the arma¬ 
ture cross field is also symmetric with respect to the point between the 
two brushes. The radial flux density of the airgap due to the field 
winding, the armature cross field and the resultant flux density are 
shown in Figure 5-12 a. As can be observed, a remaining flux density 
exists under the brushes. This flux induces a voltage in the coil 
undergoing commutation with the result that the commutation process is 
nonlinear. In order to achieve linear commutation this voltage should ~ 
be equal and opposite to the direction of the voltage induced to the 
coil due to its self inductance, thereby achieving cancellation. 
Therefore, the brushes (and the cross field) must be shifted as shown in 
Figure 5-12 b. When the magnetic vector potential is known, at a cross 
section of the motor from a finite element or other numerical method, 
the radial flux density can be calculated, since: 



(5-14) 


i M 

r 3ij> 
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In Figures 5-13 and 5-14 the magnetic vector potential and the radial 
flux density are shown for the rotor surface, Figures 5-15 and 5-16 show 
the same quantities for the stator surface. The reason the flux density 
curve is not smooth is the presence of teeth on the armature surface, 
which cause the concentration of flux in them. From these figures, the 
optimal location of the brushes can be calculated. In Figure 5-17 the 
field solution is shown for the brushes shifted 30 electrical degrees. 
In Figures 5“18 and 5-19 the magnetic vector potential values and the 
radial flux density on the rotor surface are presented for the 30° brush 
shift condition. The small amount of flux assisting the commutation can 
be observed (note arrows in Figure 5-19). 




Figure 5-12 Effect of brush position on resultant field distribution 
(a) brushes at neutral axis, (b) brushes shifted opposite to the 
rotation of the motor 
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Figure 5-17 The magnetic field for brushes 
shifted 30 electrical degrees 
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Figure 5-19 Flux density at rotor surface 
Brushes at 30 electrical degrees 
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6.0 THE STARTING OF SALIENT POLE SYNCHRONOUS 
MOTORS WITH DAMPER BARS 

The practical aspects of this problem are primarily the oscillatory 
torque and the sustained currents in the damper bars. Also of interest 
are the magnitude and waveform of phase currents and the possibility of 
the motor operating as a generator during the oscillatory cycle. The 
oscillatory torque can cause the destruction of gears and torsional 
failure of the shaft of the machine when the torque frequency coincides 
with one of the natural frequencies of the system. The currents in the 
damper bars, if sustained over a relatively long period, can increase 
the temperature of the pole face to levels which can result in mechanical 
damage. The complex waveform of the input current, and/or the operation 
of the machine as a generator can cause disturbances to the power circuit 
and possibly undesired operation of protective devices. The induced 
current in the field winding can cause elevated temperature, but the 
effect can be controlled by connecting a resistor across the field 
terminals. This resistor must be properly selected based on the prob¬ 
lems of controlling the voltage across the field terminals and current 
heating effects, and the performance of the machine must be evaluated 
with this value of resistor. 

6.1 The Electrical Equations of the Salient Pole Synchronous Machine 

Let Ui Q be the synchronous electrical speed of the machine, i.e. the 
angular frequency of the voltage applied to the terminals: 
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If the machine was not in the starting period but operating normally, 

I 

the angular velocity of the rotor would be id : 



where P is the number of poles of the machine. During the starting of 

the machine, the rotor is revolving with a speed, u> , lower than id 

’ a r * m’ mo 



where s is the slip. Defining uthe electrical speed of the rotor: 



This quantity, u) e , is fictitious and does not represent the mechanical 
angular velocity of an equivalent two pole machine operating at the same 
slip. 

Figure 6-1 depicts the three phase salient pole synchronous machine 
with a field winding, f, and damper bars. In this investigation the 
number of bars per pole has been arbitrarily chosen as 9 without signi¬ 
ficant loss in generality. 




Figure 6-1 Salient pole synchronous machine 

Kirchoffs equationsfor the machine circuits are: 

V — ip + R • X = L • I + L • I + R • I (6-5) 

where V is the vector of the voltages at the terminals of each of the 13 
current paths (3 phases, 1 field, 9 damper bars), tp is the vector of the 
flux linkages of the windings, L is the inductance matrix, R the resist¬ 
ance matrix and I.the current vector. 




100 


The entries of the inductance matrix are denoted by L.. and of the 

ij 

resistance matrix by , where i,j = a,b,c,f ,1,.. . ,9. The values of 
self and mutual inductance of the rotor circuits, ., i,j = f,l,...,9, 
are constant, due to the fact that the stator is cylindrical, and the 
self and mutual coupling does not vary with the angle 0. However, the 
self and mutual couplings between the stator circuits and between the 
stator and rotor circuits vary with the angle 0 due to the saliency of 
the rotor. These inductances can be expressed as: 

L = L + L „ cos 2 0 

aa aao aa2 

L = L + L „ cos 2 (0-120°) 

bb -aao aa2 

L = L + L _ cos 2 (0+120°) ~~ 

cc aao aa2 ■ 

L ab * - [L abo + L aa2 «*• 2 < a+30 °>] 


he ' - [L abo + L aa2 cos 2 (0 - 9O °> 1 


L ca ' - [L abo + L aa2 2 < 9+150 ">] 


L r = L .. cos 0 
af afo 


( 6 - 6 ) 


L, - = L _ cos (0-120°) 
bf afo 


L , = L , cos (0+120°) 
cf afo 


L . = [L ,. cos 0 - L . sin 0] 
ai adi aqi 

L, . = [L cos (0-120°) - L .sin (0-120°)] 
bi adi aqi 

L = [L cos (0+120°) - L .sin (0+120°)] 
ci adi aqi 


i = 1.9 


The terms L , L . are the components of stator self and mutual 
aao abo r 

inductance, and they are constant; I' aa 2 » ^afo are max; '- n,uni excursions 
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of inductances from their average values. L a ^ and I a ^ do not represent 
actual inductances; they are constants which have to be evaluated together 


with L , L , , L ~ and L , . 
aao’ abo’ aa2 afo 

The angle 0, as shown in Figure 6-1, is defined as the angle between 


the axis of the rotor and the axis of phase a in the equivalent two pole 
machine: 


0 = 0 ) t+0 = 0 ) (1-s) t + 0 

e o o o 


(6-7) 


where t is the time and 0 the initial value of 0. 

o 

From these definitions, the entries of matrix L can be calculated: 

L = -2u) L „ sin 2 0 
aa e aa2 

\t, ' - 2 “e L aa2 sln 2 (e - 120 °> 

he ' " 2 “e L aa2 sln2 < 0+12O °> 

L ab = 2u> e L aa2 sin2 ( 0+3O °) 

L bc = 2co e L aa2 sin 2 (0-9 ° O) 

L = Zo) L , sin 2 (0+150°) 
ca e aab 

L . = -a) L j. sin 0 (6-8) 

af e afo 

' '“e L afo sin <9 - 120 " ) 

L . = -a) L c sin (0+120°) 
cf e afo 


L . = -a) 
ai e 

[L adi 

sin 

0 + L . 
aqi 

cos 0] 


\i = - u e 

!L adi 

sin 

(0-120°) 

+ L . cos 
aqi 

(0-120°)] 

L . = -a) 
ci e 

[L adi 

sin 

(0+120°) 

+ L . cos 
aqi 

(0+120°)] 


The matrix R is a diagonal matrix, containing the resistivities of the. 


circuits. 
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6.2 The Effect of the Damper Bar Connections 

The damper bars in each pole are connected together at each end of 
the field pole so as to form a grill. Often this connection extends 
between poles to form a squirrel cage winding, similar to that of an 
induction motor. The type of connection has an effect on the starting 
and the transient performance of the synchronous machine, and the system 
of equations that describes it must be modified accordingly. In both 
cases the voltages of the terminals of the three phases are known and 
the voltage of the field winding terminals is zero since it is short 
circuited. The value of the field resistivity in matrix R, though, has 
to be increased by the value of the external resistor. Neglecting the 
resistance of the end ring between adjacent bars, the voltages across 
the damper bars, V^, V 2> ..., Vg have the same value, V^. 

6.2.1 Closed End Rings (Cage Winding) 

As shown in Figure 6-2, the currents and voltages at opposite poles 


are of opposite values. 



Figure 6-2 Damper winding with closed end rings 


The resistance of the pole to pole connections is taken as 2R d , yielding 
a combined resistance of R d * Thus, from Figure 6-2: 




The voltages induced in the damper bars, V. , are given by: 
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J 9c 


0 


0 


0 


• I 


0 


( 6 - 11 ) 


Since = ... = Vg = V^, equation (6-11) can be substituted for 

equation (6-5), which becomes: 



V 

c 

0 


= L 


I + L 


I + R' 


( 6 - 12 ) 



where the entries r .. of the modified matrix R a re given by 

i = a,b,c,f or j = a,b,c,f 


r' = ; r. . 

ij iJ 


and 


(6-13) 


r. . = r + R 
xj ij d 


i»j 1»•••»9 


6.2.2 Open End Rings 


The case of open end rings presents a more complicated problem. 
The determination of the voltage across the damper bars is based on the 
fact that the total current of the damper bars in each pole sums to 


zero. 
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Figure 6-3 Damper windings with' open end rings 


Neglecting again the resistivity between adjacent bars and assuming 
the same voltage V^, across all the damper bars of a pole, equation 
(6-5) can be written as: 












considerations 
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6.3.1 External Connections 

The machine is connected to the electrical network considered an 
infinite bus through an equivalent line and a power transformer. The 
inductances of these components can be calculated from the per unit 
impedance of the transformer and the short circuit apparent power at the 
end of the line. 



Figure 6-4. Equivalent circuit for the supply of the motor. 

The value of the inductance of the transmission line and the trans¬ 
former must be added to the value of L which is calculated from the 

aao 

flux maps of the machine. 

6.3.2 Self and Mutual Inductances of the Salient Pole Machine 

In order to calculate the constants of the machine inductances, the 
electromagnetic field must be solved for several excitations at various 
operating conditions. Initially the field is solved for a certain rotor 
position and the actual currents in the phase windings, the field and 
the damper bars. The solution, besides giving the magnetic vector 
potential at every point, also yields the permeabilities of all the 
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elements in magnetic portions of the machine. Then, the rotor is moved 
so that its axis coincides with the axis of phase a. By applying current 
to phase a and calculating the electromagnetic field without changing 
the permeabilities assigned to each element, the self inductance of 
phase a and its mutual inductances with the other circuits can be calcu¬ 
lated for that position of the rotor and for the saturation of the 
current paths due to the actual currents. 

The inductances thus calculated are the inductances in equations 
(6-8) corresponding to angle 0=0. 

o 

The rotor is then revolved by 90 electrical degrees and the 
previous process repeated again without changing the saturation levels. 
The same parameters are calculated again, this time corresponding to the 
angle 0= 90°. Always keeping the permeabilities of the elements the 
same as resulted from the initial field calculation, current is applied 
to the field winding and each damper bar, and the field solved for each 
case. From the information gathered from the field solutions, the 
constant values of the self and mutual inductances are then calculated 
for the specific operating condition defined by the original current 
vector. 

In addition to the external inductance, L . , the end turn induct- 

’ ext’ 

ance should be added to the calculated value of L , and the field and 

aao 

damper bar end inductances added to the corresponding inductances. 

These inductances were obtained for the case investigated here, from 

(47) 

Kilgore's formulae . It is understood that they are not accurately 
determined, but their value is too small for a significant error to be 


introduced. 
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The procedure described above was applied to a 300 HP, 6 pole 
synchronous generator with open end ring (grill). Figure 6-5 shows the 
grid used for the solution of the magnetic vector potential and Figures 
6-6 to 6-13 show the equipotential lines for currents applied to the 
field, the damper bars and phase a at zero and ninety electrical degrees. 
Due to symmetry, the field was solved for current applied to only five 
of the nine damper bars since the parameters of the rest would be the 


same. 
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Figure 6-5 The grid for a six generator 



Figure 6-7 Current in phase a, 0 
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Figure 6-14 The algorithm for the starting of a synchronous machine 
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6.. 4 The Algorithm and the Integration Scheme 

The saturation levels in the machine must be evaluated several 
times within a cycle, since, as is shown in Figures 6-23 to 6-31 the 
flux densities within the machine vary significantly during the movement 
of the rotor over one pole pitch. After the determination of the saturation 
levels, the self and mutual inductances must again be calculated as 
described in section 6.3. The system of differential equations (6.5) as 
modified in 6.3.1 or 6.3.2 to account for the voltage at the end rings 
is solved in the time domain using a predictor corrector method^^^. At 
each time step the currents are printed and prepared for plotting, and 
the integral of their square calculated. The integration is stopped 
when steady state is reached. 

Since there are no eddy currents in the iron portions of the machine, 


the torque, M, can be calculated faster and accurately by taking the partial 
derivetives of the coenergy of the field with respect to the angle: 



The flow chart of this algorithm is given in Figure 6-14. 


The currents and torque computed using this algorithm were compared 
to the results of a simulation of a start up of the same machine, based 




117 


(49) 

on test data by J.R. Misage . The phase currents and torque from 
this simulation and results of sample runs are shown in Figures 6-15 and 
6-16. 

In Figures 6-17 to 6-20, the steady state currents in the phases 

the field and the damper bars are plotted as functions of time. They 

compare favorably with the current waveforms published by Jovanovski in 
(27-30) 

1969 , from test data and analytical calculations. 

Figures 6-24 to 6-31 show the field in the machine for consecutive 
time instants when the angular velocity of the magnetomotive force is 
not the same as the angular velocity of the rotor; rather the machine is 
operating at a slip of 0.5. The time between two consecutive instants 
is one tenth of a cycle. 

Figure 6-32 shows the torque oscillations for this case and for 


other slip values. 























Figure 6~22 
The field at a cross section. Step 1 


Figure 6~23 
Step 2 



Figure 6-24 


Figure 6 - 25 
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Figure 6-32 Torque at various slips. 
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7.0 CONCLUSIONS 

The effects of saturation and eddy currents in magnetic materials 
render the conventionally defined inductances and other parameters of 
electrical machines inadequate to describe and predict the machine per¬ 
formance for modes of operation which cannot be categorized as steady 
state. The time dependent solution of the magnetic field was introduced 
in this investigation as a method for accounting for the variation, in 
time, of the machine parameters in predicting and analyzing the perform¬ 
ance of the electrical machines. 

In. order to continuously model the electromagnetic field in the 
cross section of a machine, the method of time dependent finite element 
method was used, combined with an also time dependent construction of a 
grid for the air gap region. The Maxwell stress tensor was used to 
calculate the airgap torque from the magnetic vector potential distribu¬ 
tion. Incremental inductances were defined and calculated as functions 
of time, depending on both eddy currents and saturation. The currents 
in all the machine circuits were calculated in the time domain based on 
these inductances, which were continuously updated. The method was 
applied to a chopper controlled DC series motor used for electric vehicle 
drive, and to a salient pole synchronous motor with damper bars. The 
simulation results were compared to experimentally obtained ones. 

This technique, of continuously modeling the electromagnetic field 
in electrical machines was shown to provide an insight in the operation 
of machines and account for the effect of the eddy currents and saturat¬ 
ion on the overall performance. The cases examined here do not exhaust 
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the capabilities of the method; neither do they define its limitations. 
The time dependent finite element method, although it requires repetitive 
solution of a system of linear equations, was shown to be a very powerful 
tool in taking into account nonlinearities and complex waveforms, and 
becomes more attractive when combined with a fast iterative solution 
scheme which reduces the overall computer time. 

The application and suitability of such a scheme, i.e. the pre¬ 
conditioned conjugate gradient method, was investigated. The method was 
applied to the system of linear equations resulting from the finite 
element method and it was shown that it combines both the speed required 
for the solution of time dependent problems and the ability to handle 
the ill conditioned matrices resulting from elements in iron and air 
domains. 

This work, as noted earlier, does not cover all the possible appli¬ 
cations and aspects of the continuous modeling of the electromagnetic 
field. The accuracy of the solutions for the cases already examined can 
be increased by using a finer grid, both in the solid iron regions, to 
improve the calculation of eddy currents, and in the airgap to provide a 
finer distribution of the flux densities resulting in a more accurate 
evaluation of the torque through the Maxwell stress tensor. 

The techniques used here can be applied to the skin and eddy current 
effects in solid conductors, either by constructing a finer grid for 
each of them, or by solving separately a small grid for one conductor 
bundle or one slot, every time the larger grid is solved. Application 
to separately excited or shunt DC machines and to the study of commuta¬ 
tion, would be a further extension of this work, since two or more 



circuits inductively coupled would be involved. Another possible 
application would be the case of solid salient pole machines during 
starting or transient mode. The calculation of the field, in many small 
machines in which neither the rotor is cylindrical nor the stator 
uniform, would be an interesting application of the time dependent 
computer generated grid for the airgap. 

A final improvement to the techniques described here, would be the 
application of the Newton-Raphson method to the calculation of the 
permeability of the iron. The underrelaxation method was preferred here 
merely because of the simplicity of the involved programming, but at the 
expense of computational time. 
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